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Abstract 

We investigate the compression of nuclear matter in relativistic hydrodynamics. 
Nuclear matter is described by a a — w-type model for the hadron matter phase and 
by the MIT bag model for the quark-gluon plasma, with a first order phase transition 
between both phases. In the presence of phase transitions, hydrodynamical solutions 
change qualitatively, for instance, one-dimensional stationary compression is no longer 
accomplished by a single shock but via a sequence of shock and compressional simple 
waves. We construct the analytical solution to the "slab-on-slab" collision problem over 
a range of incident velocities. The performance of numerical algorithms to solve rela- 
tivistic hydrodynamics is then investigated for this particular test case. Consequences 
for the early compressional stage in heavy-ion collisions are pointed out. 
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1 Introduction 



To investigate heavy-ion collision dynamics in realistic, i.e., 3+1-dimensional, situations by 
means of hydrodynamics requires numerical schemes to solve the hydrodynamical equations 
of motion. These schemes should reproduce analytical solutions, as far as such exist at all. 
In a previous paper we have presented two algorithms for ideal relativistic hydrodynam- 
ics, the SHASTA, a flux-corrected transport algorithm 0, and the relativistic HLLE, a 
Godunov-type algorithm ||. We have investigated their performance for the expansion of 
matter into the vacuum, where they are confronted with two problems generical for simu- 
lations of heavy-ion collisions. These are the presence of vacuum itself and the qualitative 
change in the type of the hydrodynamical expansion solution if phase transitions occur in 
the equation of state (EoS). 

As was explained in detail in JE], || , matter that undergoes a first order phase transition 
may exhibit thermodynamically anomalous behaviour in a certain range of independent 
thermodynamic variables, signalled by a change of sign of the quantity 



d 2 p 



1 - r 2 
e + p 



where c 2 = dp/de\ a is the velocity of sound, e,p, and a are the energy density, pressure, 
and specific entropy, respectively. For thermodynamically normal (TN) matter, S > 0, for 
so-called thermodynamically anomalous (TA) matter X < [[[], |J. 

As was shown in f|, for the one-dimensional expansion of TN matter a simple 
rarefaction wave is the stable hydrodynamical solution, while for TA matter such a wave 
is unstable with respect to formation of a rarefaction shock wave. Moreover, an initial 
discontinuity is bound to decay into a simple wave in TN matter, but cannot do so if 
matter is TA (or even if E only vanishes instead of becoming negative). Thus, rarefaction 
discontinuities form the stable hydrodynamical solution in the expansion of TA matter. 

On the other hand, for the compression of TA matter a compressional simple wave forms 
the stable hydrodynamic solution, but it is unstable in TN matter with respect to formation 
of a compressional shock wave. Analogously, a compressional shock discontinuity is bound 
to decay into a compressional simple wave in TA matter, while it cannot do so in TN matter. 
Thus, compression shocks are stable in the compression of TN matter. 

As outlined in |jj , in realistic cases the EoS has both TN and TA regions. Consequently, 
the hydrodynamical solution for the one-dimensional expansion is more complicated, involv- 
ing a sequence of simple waves, regions of constant flow, and discontinuities. It was shown 
in H that the same holds for the compression. The corresponding hydrodynamical solution 
was explicitly constructed in the case of a one-dimensional "slab-on-slab" collision. (For TN 
matter this test problem and the corresponding performance of the SHASTA and relativistic 
HLLE algorithms was studied in ||.) The nuclear matter equations of state used in [|J 
featured a first order phase transition between the quark-gluon plasma (QGP), described 
by the MIT bag model, and hadronic matter, described by phenomenological equations of 
state H 0. 
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In this work we first construct a nuclear matter EoS similar to that of Ref. [Q. We 
also use the MIT bag EoS || for the QGP phase, but for the hadronic phase, we employ 
a version of the a — cu-model |J (plus massive thermal pions) which features more realistic 
values for the ground state incompressibility and the effective nucleon mass than the original 
version proposed by Walecka 0. For this EoS, we construct analytically the hydrodynamical 
compression solution as described in [Q. With a tabulated version of the EoS, we then test the 
ability of the numerical algorithms to reproduce this solution. This does not only yield insight 
into the usefulness of these algorithms for applications in heavy-ion collisions and dynamical 
studies of QGP formation. In a sense it also represents an independent "numerical" proof 
for the theoretical arguments of Ref. [Q about the form of the hydrodynamical solution for 
compression of TA matter. 

This paper is organized as follows. In Section 2 we present our nuclear matter EoS. 
In Section 3 we analytically construct the hydrodynamic solution for a one-dimensional 
"collision" of semi-infinite slabs. Section 4 comprises our numerical results. In Section 5 we 
comment on one-dimensional collisions of finite nuclei and draw conclusions for numerical 
simulations of heavy-ion collisions. Section 6 concludes this work with a summary of our 
results. An Appendix contains results for modifications of the standard algorithms used 
throughout Ref. QTJ] and this work. We use natural units h — c — ks — 1- 



2 The nuclear matter EoS 
2.1 Hadron matter 

For hadronic matter we use a modification of the original a — cj-model , which was pre- 
sented in Ref. ||. For the original o — cj-model, the EoS, i.e., the pressure p as a function of 
the independent thermodynamical variables temperature T and baryochemical potential /z, 
can be rigorously derived from the a — tu-Lagrangian employing the mean-field (or Hartree, 
or one-loop) approximation of quantum many-body theory at finite temperature and density 
[T0[1 , see for instance Jll|] for an explicit derivation. A more general class of thermodynami- 



cally consistent, phenomenological equations of state for (non-strange) interacting nucleonic 
matter is defined by fll2f 



Here 



Phad(T,v) = p N (T,u;M*)+p N (T,-u;M*) + ^(^771^ (2) 

i 

rn rps 

+ nV{n)- V(n')dn'-p s S(p s ) + S(p' s )dp' s . 
Jo Jo 

p N (T, u;M*) = T^J d 3 k In [1 + exp {- (E* k - v) /T}} (3) 

is the pressure of an ideal gas of nucleons (spin-isospin degeneracy = 4) moving in 
the scalar potential S and the vector potential V R. These potentials generate an effective 

1 More accurately, V is the zeroth component of the vector potential V^, the spatial components of which 
vanish in a homogeneous, isotropic system. 
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nucleoli mass 

M — -> M* = M- S(p s ) , (4) 

where M = 938 MeV is the nucleoli mass in the vacuum, and also shift the one-particle 
energy levels 



E k = Vk 2 + M 2 — > E* k + V{n) = ^k 2 + (M*) + V(n) . 
The vector potential is conveniently absorbed in the effective chemical potential 

u = fi- V(n) , 



(5) 



(6) 



giving rise to the interpretation of as the pressure of an ideal gas of quasi-particles with 
mass M* and chemical potential v. Furthermore, 



Pi(T] mi) = -T j^- J d 3 q In 1 - exp j-^q 2 + m 2 /T 



(7) 



is the pressure of an ideal gas of mesons with degeneracy and mass m.j. In the following, 
we will only include pions (g n = 3, m, = 138 MeV), since they are the lightest and thus 
most abundant mesons. (To be consistent, one should also include thermal contributions for 
the a and to meson jTT|. However, since these mesons are much heavier than the pion, one 
can neglect their contribution for the present applications.) n is the (net) baryon density, 



n(T,fx) = — 



dp 
dfx 

_9n_ 

(2nY 



d 3 k 



1 



1 



and 



Ps(T,fjt) 



d 3 k 



M* 



e {El-u)/T + 1 
1 



e (E* k +u)/T + 1 
1 



+ 



e {E* k -v)/T + l 1 e (E*+u)/T + 1 



(9) 



is the scalar density of nucleons. It is assumed that the vector potential (which transforms 
like the zeroth component of a vector) depends only on the (net) baryon density (which 
also transforms like the zeroth component of a vector, namely the (net) baryon current iV M ), 
while the scalar potential (which is a Lorentz-scalar) depends only on the (Lorentz-) scalar 
density p s . 

Once the pressure is known, the entropy and energy density can be obtained from the 
thermodynamical relations, 



dp 
df 



e = Ts + fin — p. 



One now has to specify the potentials V, S. The choice 

V(n) = C 2 v n, S{p s ) = C 2 sPs 



(10) 
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reproduces the original a — cu-model [II]. When fitting the two parameters C v , C| to 
reproduce the experimental values for the ground state binding energy of infinite nuclear 
matter, B = 16 MeV, and the ground state density, no = 0.15891 fm -3 [|13[], the effective 
mass in the nuclear ground state is Mq ~ 0.543 M which is too small (the experimental 
value ranges around 0.7 M Jl4|), and the incompressibility is K = 9dp/dn| no ~ 553 MeV 



which is too large by about a factor of two To adjust this shortcoming of the model, 
it was suggested to introduce self-interaction terms, ~ a 3 , cr 4 , for the scalar a meson field 



in the a — cj-Lagrangian []X5||. This nonlinear a — cj-model has two additional parameters 
which enable one to also independently adjust Mq and K to the experimentally observed 
values. 

Alternatively, in Ref. M it was suggested to use 



V(n) = C v n - C 2 d n 1/3 



S(Ps) =C 2 sPs ■ 



(12) 



Although this choice introduces only one additional parameter, fitting the ground state 
incompressibility gives simultaneously reasonable values for the effective mass, for details see 
0. In the following we use the parameters C\ 
0.183, which leads to MX = 0.635 M, K = 300 MeV. 



238.08 GeV" 2 , Cf 



296.05 GeV" 2 , C\ 



2.2 Quark— Gluon Plasma 

For the QGP EoS we employ the standard MIT bag model for massless, non-interacting 
gluons and u, d quarks, i.e., 



p QGP (T, u) = — T 4 + - u 2 T 2 + -J— /i 4 - B , 

l J QLrF\ ,H>J go g f* 1627r 2 ^ 

and other thermodynamical quantities follow again from the relations 
does not depend explicitly on n for this EoS, 

P = \{e-AB) . 



(13) 

lQD. Note that p 



(14) 



For the bag constant, we use the value B = (235 MeV) 4 throughout this work. This value 
results in a phase transition temperature of T c ~ 169 MeV at vanishing baryon density (see 
below), which is roughly in agreement with what lattice QCD simulations predict |T7| . 



2.3 Gibbs Phase Equilibrium and Thermodynamical Phase Dia- 
grams 



The QGP EoS (]13D is matched to the hadronic EoS @ with (0) via Gibbs' conditions for 
(mechanical, thermal, and chemical) phase equilibrium, 



Phad — PQGP 



■ had 



TqgP , ^had — l^QGP , 



(15) 
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which leads to a phase boundary curve T*(/i*) in the T — \x plane defined by the implicit 
equation Phad(T* , /j,*) = pqgp(T* , //*) , see Fig. 1 (a). Along this curve, one can calculate the 
phase boundary values for other thermodynamical variables as a function of T*. Fig. 1 (b) 
shows the phase boundaries n H) q{T*) in the T — n/n phase diagram, Fig. 1 (c) the phase 
boundaries €h,q(T*) in the e/e — T diagram, and Fig. 1 (d) e#(7i/r(T*)) and eq{riQ{T*)) in 
the e/e — n/n diagram, respectively (e = (M — B Q )n is the ground state energy density). 
The phase transition constructed via (fL5|) is of first order, leading to a mixed phase of QGP 
and hadron matter and to a latent heat, as can be seen in Figs. 1 (b,c,d). The T-axis of 
Fig. 1 (b) maps onto the (dotted) curve e(T,n = 0) in Fig. 1 (c). Correspondingly, the 
e/eo-axis in Fig. 1 (c) maps onto the (dotted) curve e(n, T = 0) in Fig. 1 (d). This zero- 
temperature energy density is finite for finite n due to the Fermi (and, in the hadronic phase, 
interaction) energy of the fermions in the system (nucleons and quarks, respectively). This 
curve represents the minimum energy density possible for a given baryon density. 

As described in [|], the pressure as a function of energy density and baryon density, p(e, n), 
is required for solving the hydrodynamical equations numerically. It is convenient to tabulate 
this function, since calculating the pressure in the hadron phase for given e, n requires a triple 
root search (a double root search to find T, /i for given e, n, while simultaneously solving 
a fixed-point equation for the effective nucleon mass). The computational effort would be 
prohibitive if one tries to perform this "on-line" for each hydrodynamic cell in each time 
step of the hydrodynamical transport algorithm. 

Therefore, we discretize the e — n plane of Fig. 1 (d) in the range < e/eo < 20, < 
n/riQ < 12 with equidistant grid spacing Ae = to/10, An = no/20 to form a 201 x 240 mesh. 
This grid spacing proves to be sufficiently small to use the resulting table also for analytical 
calculations, see Section 4. The phase boundaries e^Q, as well as the e(n,T = 0)-curve 
are mapped onto the 240 grid points of the n/n -axis {ch,q is taken to be zero for n— values 
exceeding those where the phase boundaries meet the e(n, T = 0)-curve). 

Now values for the pressure are assigned to the mesh points. In the hadronic phase 
(defined by e(n, T = 0) < e < €h{ti) for given n) this is done using eq. (Q) and the mentioned 
triple root search, and in the QGP phase (defined by max{eg(n), e(n, T = 0)} < e < 20 eo 
for given n) eq. flUD is employed. In the mixed phase (max{e#(n), e(n, T = 0)} < e < eg(n) 
for given n) the pressure is calculated as follows. For given T*, the values of energy and 
baryon density read 



where Aq is the fraction of volume the QGP phase occupies in the mixed phase. Vice 
versa, for given e, n these equations yield values for Aq, T*. Eliminating Aq one obtains 
a single equation for T* as a function of e, n. This is solved numerically using the values 
£h,q(T*), n HiQ (T*) of Figs. 1 (b,c). Once T* is known, /x* follows from Fig. 1 (a), which 
consequently yields the pressure from Gibbs' phase equilibrium conditions flI5p . 

Similarly, every other thermodynamic quantity can be calculated as a function of e and 
n. Of particular importance are temperature, baryo-chemical potential, and entropy density. 
Once the first two are known, the entropy density is obtained from the second equation (|T0D . 



n 



e 



\ Q e Q (T*) + (l-\ Q )e H (T*) , 
\ Q n Q {T*) + (1 - \ Q )n H (T*) , 



(16) 
(17) 
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In the hadronic phase, T and \i emerge naturally from the triple root search, while in the 
mixed phase T* is obtained as described above and //* then follows from Fig. 1 (a). In the 
QGP, we distinguish two cases. For n = we also have fi = 0, and then we obtain from 
(|IJ, [n§ the simple formula T = [30(e - J B)/37vr 2 ] 1 / 4 . For finite ji we eliminate T from the 



equation of e using the equation for n. This results in a sixth order equation in fx, 

A A QQQ7T 2 

which has to be solved numerically. After that, T = [9n/2/i — fi 2 /9tc 2 ] 1 ^ 2 (which follows from 
the equation for n). 

In Fig. 2 we show (a) pressure, (b) temperature, (c) baryo-chemical potential, and (d) 
entropy density as functions of e and n, as they emerge from the above described procedure. 
Note in Fig. 2 (a) that on account of (|14]) the pressure is independent of n in the QGP phase. 
Also, in the mixed phase for small n the pressure is (almost) constant as a function of either 
e or n. The thermo dynamical identity 



2 d P 



dp 



n dp 
e + p dn 



(19) 



then implies that the velocity of sound is very small in this region. For n = 0, the pressure 
(as well as the temperature, see Fig. 2 (b)) is constant in the mixed phase, and the velocity 
of sound vanishes identically (cf. the EoS for TA matter in |l|). 

In contrast to the diagrams for the thermodynamical variables p, T, fi which enter Gibbs' 
phase equilibrium conditions (|15D, there is no obvious sign for the mixed phase (such as 
edges due to the phase boundaries) in the entropy density. The reason is that this quantity 
is a linear interpolation between sh and sq (similar to (|16|, [IT])) in this phase. 

For the hydrodynamical simulations intermediate values of thermodynamic quantities are 
calculated by two-dimensional linear interpolation. If e > max{20 eo, e(ri, T = 0)} for given 
n, the QGP EoS (0) is used directly to calculate the pressure, and temperature, chemical 
potential, and entropy density are obtained as described above. 



3 Compression of nuclear matter 

In this section we construct the analytical solution to the one-dimensional compression of 
nuclear matter with the EoS presented in Section 2. The initial condition of the "slab-on- 
slab" collision problem is 

e(x, 0) = e , — oo < x < oo , (20) 
n(x, 0) = n , — oo < £ < oo , (21) 

v(x,0) = -°°<*< > (22) 

v ' I -v C m , < X < oo . v ' 
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In a real heavy-ion collision (see also Section 5), vcm corresponds to the velocity of the 
nuclei in the "equal-velocity frame". For identical nuclei, this is equivalent to the center- 
of-momentum (CM) frame. In the subsequent time evolution, momentum conservation will 
force matter to stop at x = and energy conservation to convert its kinetic energy into 
internal energy, leading to compression waves travelling symmetrically away from the origin, 
leaving in their wake compressed and heated matter at rest. For small vqm, the energy 
and baryon number densities obtained in the compressed state correspond to hadronic, i.e., 
TN matter. Thus, compression proceeds via single shock waves. Energy, momentum, and 
(net) baryon number conservation across the (space-like) shock discontinuity read in the rest 
frame of the shock (quantities with subscript "0" correspond to the uncompressed state) 

T oi = T oi ^ T u = T u ^ N i = N i ^ ^ 

where = (e + p)u^u u — pg^ u is the energy momentum tensor for an ideal fluid, and 
= nu^ is the (net) baryon number current. The final compressed state (e, n) obeys the 
Taub equation [|18| 

wX - w X - (p - po) {X + X ) = , (24) 

where w = e + p is the enthalpy density and X = w/n 2 the generalized volume^. The Taub 
equation defines the so-called Taub adiabat in the p—X plane. It is the set of all possible final 
states which are in accordance with energy, momentum, and baryon number conservation 
across a single shock discontinuity. A particular final state is fixed by specifying the involved 
velocities. For instance, in our case the velocity vcm of uncompressed matter in the rest 
frame of compressed matter is given. On the other hand, in the rest frame of the shock, 
uncompressed matter flows into the shock with velocity Vq and compressed matter emerges 
with velocity v. These velocities obey |TJ| 



(p-PoXe + Po) v 2 _ (p-Po)(eo+p) ^ 



" (e-e )(e +p) ' (e - e )(e + p ) 

Obviously, vcm = \vq — v\/(l — vqv) and thus 



2 (e + p)(e + p ) 
1cm — 

WWq 



(e + po) n 
Wq n 



e/n 



(26) 



since po = in the ground state of nuclear matter. In deriving this equation we have made 
use of the Taub equation (p4|) . Given vcm , this equation fixes e/n, or with the Taub equation, 
e and n in the final compressed state. 

Fig. 3 (a) shows the Taub adiabat with the center (X ,p ) = (e o /no,0) for the EoS of 
Section 2. The different parts belonging to final states in the hadron, mixed, and QGP phase 
are marked correspondingly. One observes a region of final states (A— C) where the chord 
connecting them with the center intersects the Taub adiabat. For these final states a single 



compression shock is not the hydrodynamically stable solution 0, |12|, [19|]. 



2 Note that this definition differs from the one given in pj] for the (net) baryon-free case. 
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To see this and to determine the stable solution one first notes that point A is a Chapman- 
Jouguet (CJ) point of the Taub adiabat. At this point, the specific entropy a = s/n assumes 
a (local) maximum and the slope of the Taub adiabat and the Poisson adiabat (the curve of 
constant specific entropy) are identical |L9] 



dx 

dp 



dx 

dp 



(27) 



TA 



Furthermore, above the CJ point the chord connecting a final state {X,p) with the center 
has a smaller slope than the Taub adiabat at (X,p), i.e., 



dX 
dp 



< 



TA 



x-x 

P-Po 



P > Pcj 



(28) 



with the equality holding at the CJ point. Subtracting this equality from (]28|) , we obtain 

dX/dp\ TAiP - dX/dp\ a>Pcj 



P-Pcj 

In the limit p — ► pcj we conclude using fl2"7|): 

d 2 x 



< 



P> Pcj ■ 



(29) 



dp 2 



= £ < for p > pcj 



(30) 



<T,PCJ 



From thermodynamical identities, S = H/n 2 c®, and therefore, above the C J point matter is 
TA. This implies that, once the CJ point is reached via a single compression shock, further 
compression can be achieved only through a compressional simple wave. The hydrodynamical 
solution will therefore consist of a single compression shock to the CJ point and, attached 
to it, a compressional simple wave (see Fig. 5 below). Since matter emerges from the shock 
with the velocity of sound (since the compressed state corresponds to the CJ point |T9"|), and 
since for a simple wave the matter velocity relative to the wave profile is also the velocity 
of sound, the simple wave profile does not move relative to the shock, and they will remain 
attached to each other in the stationary state. Since the specific entropy is constant for 
continuous solutions of ideal hydrodynamics, final states of compressed matter lie on the 
Poisson adiabat through the CJ point, rather than on the Taub adiabat, see Fig. 3 (b). 

Compression of matter can proceed through this configuration only until the point of 
inflection, d 2 X/dp 2 \ a = 0, on the Poisson adiabat is reached (point B in Fig. 3 (b); this is 
usually the point where the Poisson adiabat leaves the mixed phase and enters the QGP 
phase). Since matter becomes TN at this point, further compression is achieved by another 
compression shock attached to the head of the compressional simple wave (see Fig. 6 below). 
This shock is also described by the Taub equation (f2~4|), however, X , p , and w have to 
replaced by the corresponding quantities X, p, and w at the head of the compressional simple 
wave. The final compressed state lies in the pure QGP phase and is determined as follows. 
For a stationary configuration, the (baryon) flux iV 1 = rvyv at the head of the compressional 
simple wave has to be identical to that through the shock front. In the rest frame of the 
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latter it is determined by (N 1 ) 2 = — (p — p)/ (X — X) (which follows from the conservation 
laws (|23|) ; the interpretation is that the square of the baryon flux through the shock equals 
the negative slope of the chord connecting initial and final state of the shock). On the other 
hand, the matter velocity at a fixed point of a simple wave profile is equal to the velocity 
of sound, i.e., in the rest frame of that point the baryon flux is iV 1 = nc s /(l — c 2 ) 1 ^ 2 . Now 
consider this point to be the head of the simple wave, or equivalently, the position of the 
compression shock, if the configuration is stationary. Since dp/dX\ a = —n 2 c 2 /(l — c 2 ) (which 
follows from thermodynamical identities), we arrive at the condition 



dp 
dX 



a.X X ~ X 



P ' P (31) 



for a stationary hydrodynamical configuration. This means that the final compressed state 
(X,p) on the Taub adiabat with center (X,p) is determined by the intersection of this 
adiabat with the tangent to the Poisson adiabat at (X,p). 

In order to reach higher compression through this shock, one has to decrease the slope 
i.e., the state (X,p) has to "move backwards" along the Poisson adiabat towards the CJ 
point A. The compressional wave will consequently become "shorter" and the compression 
shock stronger. This proceeds until the CJ point A is reached. Then, the compressional 
simple wave vanishes and the shock (O-A) merges with the shock (A-C). Since the chord 
(A-C) has the same slope as the chord (O-A), the baryon flux through the two shocks 
is identical and therefore, the single shock (O-C) becomes the hydrodynamically stable, 
stationary configuration. Above C, the single shock described by the Taub adiabat of Fig. 3 
(a) represents again the stable hydrodynamical solution. 

In practice, one obtains part (B-C), the so-called wave adiabat H, when calculating 



the Poisson adiabat: for each (X,p) one simply solves (31) simultaneously with the Taub 
equation with center (X,p). The set of all final compressed states (X,p) thus obtained was 
called generalized shock adiabat in Ref. |4[]. 

Let us now complete the solution of the hydrodynamical problem. As was the case for 
the expansion of semi-infinite matter studied in the stationary hydrodynamical solution 
is of similarity type, i.e., its profile is constant as a function of ( = x/t. Due to causality 
and the symmetry of the problem, it suffices to consider the range < ( < 1. Below and 
at the CJ point, compression proceeds through a single shock, travelling to the right with 
velocity v s h = —v, where v is given by fl25|). Given vcm, the final state energy density e and 
baryon number density n are obtained as solutions of (f24|) and (|26|) . Other hydrodynamical 
variables can be inferred from the EoS and the definition of T^ u and N 11 . Profiles of T 00 /eo 
and T as functions of ( are shown in Fig. 4 for vcm = 0.7. 

The CJ point itself has to be determined numerically by the requirement of maximum 
specific entropy. For the above EoS, p C j — 1.941 e , X c ,j ^ 0.422 X , (e c ,j — 6.971 e , n C j ^ 
4.596 n ,) with a specific entropy a ~ 3.545. For further purpose, let us denote the value of 
vcm required to reach the CJ point by vcj- From ( [26]) we infer vcj — 0.752. 

Above the CJ point and below point B, compression proceeds through a shock to the 
CJ point and an attached compressional simple wave. The shock moves with v s h = [v' sh — 
u cm]/[1 — v 'sh v cm] where v' sh is the shock velocity in the rest frame of incoming matter. The 
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latter is identical with — vq, cf. (BH) , i.e., v' sh ~ 0.878. Obviously, matter emerges from the 



shock with velocity vcj in the rest frame of incoming matter, i.e., in the CM frame matter 
emerges from the shock with velocity [vcj — «cm]/[1 — vcjVcm]- This is, of course, still 
negative above the CJ point, since the subsequent compressional simple wave will further 
decelerate matter. 

Consider the compressional simple wave moving to the left (i.e., that at x > 0) in the 
rest frame of matter emerging from the shock. Constancy of the Riemann invariant 7£_ for 
this wave JTJ implies that the velocity as a function of the energy density reads 

< M ^anh{£ dE '^}. (32, 

The integral is calculated via a centered Riemann sum. The required values of p(e') and e' 
are taken from Fig. 3 (b). In that calculation it was also necessary to locally determine the 
tangent to the Poisson adiabat (in order to construct the wave adiabat (B-C)). The slope 
of the tangent is dp/dX\ a = — n 2 c 2 s /{l — c 2 s ) and this relation can be easily inverted to yield 
the local velocity of sound c s (e') required for the integration in fl32|). 

In the rest frame of incoming matter, the velocity of matter on the simple wave is 
v ' w ( e w) = [v'l(e w ) + v C j]/[l + v^(e w )v C j], and consequently in the CM frame v w (e w ) = 
[ v ' w ( e w) ~ v cm}/[^ — v ' w ( e w ) v cm\- The complete simple wave profile is obtained by succes- 
sively increasing e w in ([32] ) until v w becomes zero. Then, the simple wave has completely 
decelerated matter in the central region. 

The similarity variable ( corresponding to a given e w is determined via |I| 

£ Vwjtw) ~\~ c s\£w) (33) 

1 + v w (e w ) c s (e w ) 

The final energy density e where matter comes to rest determines the position £ = c s (e) of 
the head of the simple wave. Profiles of T 00 /eo and T as functions of ( are shown in Fig. 5 
for vcm = 0.8. 

Point B on the generalized shock adiabat is numerically determined as ps — 2.478 eo, 
Xb — 0.198Xo (cb — 18.271 eo, Ub — 10.236 no). In the rest frame of matter emerging from 
the shock (O-A), the velocity at the head of the simple wave is ^(es) as given by (|32"D, and 
thus, in the rest frame of incoming matter, vb = + u cj)]/[1 + v w( e B) vcj] — 0.820. 

Therefore, the CM velocity required to reach point B is vqm = vb- 

Above point B and below C, the compression profile consists of a shock to the CJ point 
A, a compressional simple wave, and a second shock. The position of the first shock and 
the respective hydrodynamic quantities are determined as above. In the rest frame of the 
second shock, the velocities obey the relations (p5|), except that quantities with subscript 
have to be replaced by the quantities e, p, and v. In the rest frame of compressed matter, 
the second shock moves with velocity 



Vsh = -V 



(P-P)(P + £) 1 V2 

(e -e)(p + e) ' [6V 
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while matter flows into that shock with velocity 



v,„(e = 



vv 



(p-p)(e-e) 
(p + e)(p + e) 



1/2 



(35) 



Using ( |3"2"D (after boosting it to the global CM frame) for the left hand side, this equation 
represents a condition to determine e, p. In practice, one constructs the compressional simple 
wave starting from the CJ point as in the previous case, but simultaneously checks for each 
v w {^w) whether the condition (|3~5p is fulfilled. At first, the left hand side will be larger than 
the right due to the following reason: a compressional wave which is too "short" implies a 
shock that is too strong (the corresponding final state lies above the true final state on the 
wave adiabat), and a stronger shock is more effective in decelerating matter. For a given 
vcmi a shock that is too strong will result in a net positive velocity of matter in the final 
compressed state, instead of bringing it only to rest, or in other words, v w {e) was not yet 
sufficiently negative in order that matter is just stopped by the second shock. Analytic 
profiles for the hydrodynamic variables are shown in Fig. 6 for vqm — 0.825. 

Point C is numerically determined aspc = 2.703 e , X c = 0.195 X (ec = 18.944 e , n c = 
10.525 no). From this point onwards, single shock solutions are again hydrodynamically 
stable. Thus, from fl26|) we determine the CM velocity to reach that point as vqm = v c — 
0.831. Profiles for the hydrodynamic variables are presented in Fig. 7 for vcm = 0.9. 

This concludes the discussion of the analytic construction of the compression wave pro- 
files. We finally note that in case point A of the Taub (or generalized shock) adiabat is not a 
true CJ point but a kink in the adiabat^, also a double shock solution is possible for certain 
values of vcm [f§- Since this situation does not occur for the EoS considered here, we do not 
discuss it in detail. 



4 Numerical results 

In this section we present numerical solutions for the compression of nuclear matter. The 
first test is whether the tables of thermodynamic quantities constructed as described in 
Section 2 are sufficiently accurate to reproduce the analytical solution presented in the 
preceding section. To this end, we solved the Taub equation ( p4|) and constructed the 
generalized shock adiabat using linear interpolation on these tables. The curves of Fig. 3 
are satisfactorily reproduced, even the wave adiabat, the construction of which involves the 



simultaneous solution of eqs. ( pl|) and ([Jl]). The position of the CJ point can be found within 
1% accuracy. 

Then we checked whether the analytic shape of a simple compression wave can be repro- 
duced. To this end, a discretized form of the thermodynamical identity (|19D using tabulated 
values for the pressure was employed in the integrand of eq. (|3"2"p. Also in this case, agree- 
ment was found to be rather good, confirming that the tabulation of the thermodynamic 
quantities is sufficiently accurate. 

3 This could happen for particularly "stiff" hadron matter equations of state. Then, this point corresponds 
to the phase boundary between hadron and mixed phase matter. 
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The SHASTA and relativistic HLLE algorithms were described in detail in Ref. JTJ. To 
test their numerical performance in the "slab-on-slab" collision problem, we choose the four 
different CM velocities vcm = 0.7, 0.8, 0.825, and 0.9, for which the analytical profiles were 
explicitly constructed as discussed in Section 3. The numerical profiles of the CM frame 
energy density T 00 (normalized to the ground state energy density e ) and the temperature 
T (in MeV) are shown in Figs. 4 (a-c) for the relativistic HLLE (At/ Ax = A = 0.99) and in 
Figs. 4 (d-f) for the SHASTA (A = 0.4) for vcm = 0.7. (The presentation of quantities in 
terms of the similarity variable ( is advantageous since in this way one can easily monitor the 
approach of the numerical to the analytical solution Jl|.) The resolution of the shock front 
is rather good (~ 4 grid cells) for both algorithms already after a few time steps (~ 50). 
The SHASTA produces a small overshoot at the shock front. 

In Fig. 5 we show the corresponding results for vqm — 0.8. The approach to the an- 
alytical solution is slow (about 500 time steps for both algorithms). At early times the 
compression configuration rather resembles a shock front that is broadened by viscosity [^(J . 
The SHASTA run features a non-propagating instability on the shock plateau at x ~ 25 (as 
a function of ( = x/t, such an instability moves of course to the left in the course of time, 
cf. Fig. 5 (e)). This instability can be removed by decreasing the antidiffusion fluxes, see 
Appendix. The HLLE run shows a slowly decaying overshoot in the energy density around 
x = ( = 0. It is a remnant of the first few time steps when numerical transients cause an 
overestimate of the shock plateau (this effect is clearly visible in Figs. 4 and 7). Apart from 
these phenomena, both algorithms reproduce the position and shape of the compressional 
wave configuration quite well after about 500 to 1000 time steps. 

The origin of the instability (Fig. 5 (d)) and the overshoot (Fig. 5 (a)) and the reason for 
their persistence are easily understood considering the fact that the final state corresponds 
to TA mixed phase matter where pressure gradients are small (cf. Fig. 2 (a); the same 
holds for temperature gradients, cf. Fig. 2 (b), therefore, corresponding phenomena do not 
occur in the temperature profiles, Figs. 5 (c,f)). On one hand, this facilitates compression 
of matter (one has to exert less force in the compression) and subsequently the creation of 
local density maxima. On the other hand, the usual driving force for the expansion of local 
density maxima is absent. They can only decay on account of numerical diffusion. Therefore, 
the instability produced by the standard version of the SHASTA (Figs. 5 (d,e)) is removed 
when the antidiffusion is decreased, since this increases the numerical diffusion. On the other 
hand, the numerical diffusion of the relativistic HLLE is just large enough to slowly damp 
out the overshoot at x = ( = (Figs. 5 (a,b)). We finally note that in the mixed phase a 
density increase (along the simple compression wave) corresponds to a temperature decrease, 
cf. Fig. 5 (b,c,e,f), reflecting once more the TA nature of this phase. 

Fig. 6 shows our results for vqm = 0.825. While the relativistic HLLE is well able 
to reproduce the rather complex compressional wave configuration (modulo unavoidable 
numerical dissipation), the SHASTA fails completely: apart from a slowly growing, non- 
propagating instability at x ~ 10, a single shock front forms instead of the configuration 
consisting of two shocks and a simple wave. The reason is again that the numerical dissipation 
is too small (see also Appendix). Reducing the antzdiffusion in the SHASTA therefore 
considerably improves the reproduction of the analytical profile, see Appendix for details. 
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Finally, in Fig. 7 the profiles are shown for vcm — 0.9. Since the hydrodynamically 
stable solution is again a single compression shock, this case is similar to Fig. 4, i.e., both 
algorithms reproduce the analytical solution rather well after a few time steps. Comparing 
Figs. 7 (b) and (e), the resolution of the shock front in the SHASTA appears to be worse 
than for the HLLE. However, one has to remember that, due to the smaller CFL-number A 
for the SHASTA run, the profiles extend over fewer cells after the same number of numerical 
time steps |]. Close inspection reveals that the shock front is resolved over ~ 3 grid cells 
for both algorithms after about 30 time steps. 

In conclusion, despite the simplicity of the numerical algorithms the (partly rather com- 
plex) analytical solutions are reproduced remarkably well with the relativistic HLLE and 
also with the SHASTA after decreasing the antidiffusion. This constitutes an independent 
"numerical" proof of the correctness of the arguments presented in |4| and Section 3. It is 
remarkable that an usually unwanted feature like numerical diffusion is now necessary to 
reproduce the physical solution. The reason is that in TA matter this solution reacts sensi- 
tively to numerical instabilities. In order to suppress the latter and consequently make the 
transport algorithm more robust, one has to increase the numerical diffusion. In the case of 
single compression shocks the approach to the analytical solution requires about 50 numeri- 
cal time steps. However, to reproduce the more complex compressional wave configurations 
for TA matter takes about a factor 10 to 20 longer. Consequences for the simulation of 
heavy-ion collisions are pointed out in the next section, where the compression of finite sys- 
tems is studied. We note that for the HLLE runs we have used a constant velocity of sound 
c\ = 1/3 in the signal velocity estimates, cf. also [JJ. Using the physical sound velocity has 
already been shown to produce unphysical solutions in the expansion of matter into vacuum 
]1|]. As demonstrated in the Appendix, such problems occur also in the compression 

5 Compression of finite systems 

We consider a one-dimensional collision of two (equal) nuclei with rest frame radius R. The 
nuclei have the CM velocity vcm and are initialized in the moment of first contact. The 
initial condition for this problem is 



4 Of course, the choice c 2 s = 1/3 is only safe as along as the physical velocity of sound does not exceed 
this value. For instance, for pure hadronic matter described by the a — w-model the velocity of sound 
even approaches the causal limit for large e, n. Fortunately, for the EoS constructed in Section 2, the phase 
transition to the QGP phase (with c 2 s = 1/3) prevents this, and no problems were encountered for the test 
cases studied here. 




(36) 



(37) 



(38) 
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The compression ends when the compression front has traversed the nucleus. In the CM 
frame, this corresponds to a time J|] 

tF ^OM{l-v' sh V CM ) R 
V'sh 

where v' sh is the velocity of the compression front in the rest frame of the incoming nucleus 
as calculated in Section 3, i.e., it is either the velocity of the single compression shock in the 
case vcm < vcj or vcm > vc, or that of the shock to the CJ point for vcj < vcm < ve- 
in the latter case, v' sh ~ 0.878, see Section 3. In Fig. 8 tp (in units of R) is shown as a 
function of vqm- If vcm is less than the velocity of sound in the ground state of nuclear 
matter c Sj0 = (-^o/9/io) 1 ^ 2 (— 0.190 for the parameters of Section 2), no shock discontinuity 
can form ||19|| . The point where the generalized shock adiabat enters the mixed phase and 
also point C on this adiabat (cf. Fig. 3 (b)) are visible as kinks in tp. 

In order for the numerical algorithms to approach the correct analytical solution for times 
t < tp, i.e., before expansion sets in, the grid spacing has to be sufficiently small. We have 
seen in the previous section that in the case of single compression shocks the shock profile 
is accurately reproduced after about 50 time steps. Therefore, the grid spacing should be 
chosen smaller than Ax ~ £f/50A. To give an example, for tp ~ R, Ax < 0.02 R for the 
HLLE (A = 0.99) and Ax < 0.05 R for the SHASTA (A = 0.4). For a typical nuclear radius 
of R ~ 5 fm, a feasible grid spacing for the HLLE would therefore be Ax < 0.1 fm and for 
the SHASTA Ax < 0.25 fm. However, for ultrarelativistic collisions where tp ~ 0.2 R, Ax 
should at least be a factor of 5 smaller. 

We do not show explicit calculations for the case of compression via single shock dis- 
continuities, since the result is obvious: after the incident nuclei are completely stopped by 
the compressional shock waves, expansion sets in. The only remaining question is whether 
the performance of the algorithms for this expansion of baryon-rich nuclear matter with a 
realistic nuclear matter EoS is of similar quality as for the (net) baryon-free case studied in 
Ref. |l|. This investigation is out of the scope of the present work. 

For the complex compressional wave patterns of Figs. 5, 6 we have to good approximation 
tp ~ 0.57 R for all cases, cf. Fig. 8. We have seen that of the order of 500 time steps are 
necessary to approach the analytical profiles. Thus, the grid spacing should be chosen 
smaller than 0.001 R ~ 0.005 fm for the HLLE and smaller than 0.003 R ~ 0.015 fm for 
the SHASTA. The numerical effort to calculate a collision on such a fine grid is definitely 
prohibitive, especially for multi-dimensional problems. It is therefore of interest whether 
the distortion caused by a too coarse grid spacing really changes physical observables, such 
as particle spectra, in the final state. 

The calculation of these spectra requires assumptions about the particle emission process 
from the fluid and, ultimately, about the "freeze-out" of the fluid. This is out of the scope 
of the present work and will be subject of a forthcoming paper For the moment, 



we only study the spatial CM energy density distribution of the fluid, which serves as 
input for a "freeze-out" calculation. We compare the time evolution of this quantity for 
one-dimensional collisions of finite nuclei at vqm = 0.8, calculated with the HLLE for 
Ax = 0.01 R and Ax = 0.001 R (Fig. 9), and calculated with the SHASTA for Ax = 0.025 R 
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and Ax = 0.0025 R (Fig. 10). For the SHASTA runs we have reduced the antidiffusion 
fluxes, since this stabilizes the algorithm and improves results considerably, cf. Appendix. 
Let us note that energy, momentum, and baryon number conservation for all runs is better 
than 10~ 3 . 

In Fig. 9 one observes that for the compressional stage, Figs. 9 (a,b), the finer grid spacing 
leads, as expected, to a much more accurate description of the analytical profile, cf. Fig. 5. 
However, the final CM energy density profiles emerging after the expansion stage, Figs. 9 
(c,d), are quite similar. Indeed, the run with the coarser grid spacing even yields a smoother 
profile and does not show "terraces" as can be seen in (b,d) on the low density tails in 
either compression or expansion stage. Furthermore, small-scale numerical instabilities are 
observed around the stationary point of the profile due to the exceedingly long calculation 
time in (d). Consequently, for the HLLE the use of a grid spacing which is sufficiently fine 
to reproduce the compressional stage yields no obvious advantage when one is interested 
in final state fluid quantities, but it is even disadvantageous from the point of view of the 
required calculational effort (the same physical time corresponds to 10 times the number of 
numerical time steps, the smaller grid spacing requires 10 times more grid cells on the same 
physical length scale). 

From the physical point of view we note that the expansion of the system proceeds similar 
as in the (net) baryon-free case studied in |l[]: the final compressed state consists of TA mixed 
phase matter, and is therefore consumed by a rarefaction discontinuity. In Figs. 9 (c,d) one 
clearly observes this shock wave travelling into the high density zone. Hadronic matter is 
expelled from the shock and subsequently expands via a simple rarefaction wave. 

For the SHASTA (Fig. 10) the distorting effect of a too coarse grid spacing is more 
obvious. Clearly, not only is the reproduction of the compression profile insufficient, also the 
rarefaction shock wave in the expansion stage is no longer discernible. On the other hand, 
the run with the finer grid spacing Ax = 0.0025 R produces excellent results, although there 
exists again a tendency to produce "terraces" on the tail of the hadronic expansion wave, as 
for the HLLE run with As = 0.001 R. 

However, as far as the final state profiles and particle spectra at freeze-out are concerned, 
all runs would give similar results. It is therefore questionable whether high precision cal- 
culations with a small grid spacing are really compulsory for the simulation of heavy-ion 
collisions. We note that results for the case vqm = 0.825 are rather similar, wherefore we do 
not show them explicitly. 

6 Conclusions 

In this paper we have studied the compression of nuclear matter in one-dimensional hy- 
drodynamical "slab-on-slab" collisions. First, a nuclear matter EoS was constructed, which 
describes the hadron matter phase in terms of an improved version of the original o — lo- 
model. For the QGP phase, the MIT bag model EoS was employed, and both phases were 
matched via Gibbs' phase equilibrium conditions for a first order phase transition. In nu- 
merical algorithms to solve relativistic hydrodynamics, this EoS must be used in tabulated 
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form, since an "on-line" determination of the pressure (which is required in the solution of 
the hydrodynamical equations at various intermediate steps, cf. [0) would be calculation- 
ally prohibitive. However, the usefulness of the tables for thermodynamic quantities is not 
restricted to the hydrodynamical problems studied here, they can also be applied to future 
and, most important, more realistic multi-dimensional calculations. 

Subsequently we have constructed the analytical solution to the one-dimensional com- 
pression of nuclear matter for various incident velocities. Since nuclear matter with a first 
order phase transition is TA in certain regions of the independent thermodynamical variables, 
the hydrodynamically stable compression solution is no longer a single shock discontinuity, 
as for TN matter, but consists in general of a sequence of shocks and compressional simple 
waves. We remark that the constancy of the specific entropy on compressional simple waves 
leads to a plateau in the respective excitation function which might in turn lead to a cor- 
responding plateau in the excitation function of the pion multiplicity |22j. This would be a 
clear signal for the phase transition between hadron matter and the QGP. 

Then we have investigated the performance of numerical algorithms for this particular 
test problem. Both the relativistic HLLE and the SHASTA algorithm are well able to 
reproduce the analytical profiles, provided the numerical diffusion is sufficiently large. In a 
sense, this constitutes an independent numerical "proof" for the correctness of the theoretical 
arguments presented in Ref. [|J and in Section 3 of this work concerning the hydrodynamically 
stable compression solution. The algorithms require about an order of magnitude longer to 
approach the more complicated compressional solution in the region where nuclear matter is 
TA than they need for the reproduction of the compressional shock waves in TN matter. To 
account for this, one has to choose a sufficiently fine grid spacing Ax in collisions of finite 
nuclei at the respective incident energies, in order that the analytical profile is reproduced 
before the compression wave has traversed the nucleus and expansion sets in. However, we 
have found that the final profiles for a calculation employing a too coarse Ax do not differ 
significantly from those for a finer Ax, so that the computational effort can be considerably 
reduced, at least as long as one is interested in final state observables onlyQ. 

We anticipate that the algorithms presented in detail in Ref. ||] and studied for various 
non-trivial test problems in and the present work will find numerous applications in 
multi-dimensional hydrodynamical problems. Of special interest is a deeper understanding 
of the flow as discovered in BEVALAC experiments a decade ago |24| and just recently 



confirmed by quantitatively excellent data from the EOS-collaboration [EH]. Flow was also 



observed in recent AGS experiments pq| . At present it is unclear whether its magnitude can 
be completely accounted for in the framework of microscopic cascade models [E7] or whether 



it is consistent only with a phase transition in the nuclear matter EoS [28 



5 However, the effect wilf be negligibfe even on observables that are most sensitive to the initial hot stage, 
like photons or dileptons [p3| , since the average temperatures are also very similar for the various runs of 
Figs. 9, 10. 
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A Appendix 

In this Appendix we study the effect of modifications of the algorithms as used in the main 
part of this work. First, we demonstrate that the performance of the SHASTA in the 
reproduction of the complex compressional wave configurations of Figs. 5, 6 is considerably 
improved by decreasing the antidiffusion fluxes. As discussed in the Appendix of [p]], to this 
end the first factor 1/8 in eq. (17) of is replaced by 1/10. The effect is shown in Fig. 
11. One observes that the instabilities in Figs. 5, 6 (a,b,d,e) are removed due to the larger 
numerical dissipation introduced by reducing the antidiffusion fluxes. Contrasting Fig. 6, 
the wave configuration for vcm = 0.825 is now in reasonable agreement with the analytical 
profile. 

Finally, we discuss the effect of using the physical velocity of sound in the signal velocity 
estimates for the relativistic HLLE. In Fig. 12 we confront the correct solution for (a) vcm = 
0.8 and (b) vqm = 0.825 obtained with a constant velocity of sound c\ = 1/3 in the signal 
velocity estimates (cf. eqs. (29, 30) of Ref. |]) with the corresponding results (c,d) using 
the physical sound velocity calculated according to ([T^) on the table of pressure values. 
Comparing (a) and (c) one observes that instead of reproducing the simple compression 
wave the code now grossly overestimates the strength of the first shock which even accelerates 
matter to positive velocities. Subsequently, a second rarefaction shock brings matter to rest. 
The final state has a smaller energy density than the correct solution. 

A comparison of (b) and (d) reveals that the code fails completely to reproduce the 
configuration consisting of a compressional simple wave between two shocks. Instead, a 
single compressional shock is formed which does not decay, as it should for TA matter. The 
situation is rather similar to what the standard version of the SHASTA yields for this case 
(Fig. 6 (e)). A too small diffusion was then identified as the reason for the failure to produce 
the analytical result. This is a natural explanation also for understanding Fig. 11 (d): the 
physical sound velocity in the mixed phase is small, and a small sound velocity decreases 
the signal velocity estimates [l], and consequently, the numerical diffusion. 
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Figure Captions: 

Fig. 1: Nuclear matter phase diagram in the (a) T — fi, (b) T — n/n , (c) e/e — T, and (d) 
e/e — n/riQ plane. The dotted line in (c) corresponds to the n = 0-axis in (b), the dotted 
line in (d) to the T = 0-axis in (c). 

Fig. 2: (a) p/e , (b) T, (c) /i, and (d) s/n as functions of e/e and n/n . (For pre- 
sentation purposes, only every second meshpoint on the 201 x 240 mesh is plotted.) 

Fig. 3: (a) Single shock Taub adiabat with the center (Xq,pq) (point O) and the CJ point A. 
(b) Generalized shock adiabat. The part (O-A) and that above C are identical to the Taub 
adiabat in (a). Part (A-B) is a Poisson adiabat (B is the inflection point of this adiabat), 
part (B-C) the wave adiabat (for details concerning its construction see text). The dotted 
line is the unstable part of the Taub adiabat of Fig. 3 (a). 

Fig. 4: "Slab-on-slab" collision for vqm = 0.7 as calculated with (a-c) the relativistic 
HLLE algorithm (Ax = 1, A = 0.99) and (d-f) the SHASTA (Ax = 1, A = 0.4). (a,d) CM 
frame energy density T 00 (normalized to the ground state energy density e ) as a function of 
x for 10,20,30, 100 time steps, (b,e) T 00 /eo and (d,f) temperature T as functions of the 
similarity variable £ = x/t, for 10 (stars), 20 (squares), 30 (diamonds), 50 (dotted line), and 
100 time steps (dashed line) in comparison to the analytical result (full line). 

Fig. 5: The same as in Fig. 4 for v C m = 0.8. In (b,c,e,f) stars correspond to 50, squares to 
100, diamonds to 200, the dotted line to 500, and the dashed line to 1000 time steps. 

Fig. 6: The same as in Fig. 5 for vqm = 0.825. 

Fig. 7: The same as in Fig. 4 for vcm = 0.9. 

Fig. 8: The CM time tp the compression front needs to reach the edge of the incoming 
nucleus as a function of vqm- 

Fig. 9: The CM frame energy density (in units of e ) as a function of x (in units of 
the nuclear radius R in its rest frame) for a one-dimensional collision of finite nuclei at 
Vcm — 0.8, calculated with the relativistic HLLE (A = 0.99) for (a,c) Ax = 0.01 R and (b,d) 
Ax = 0.001 R. (a,b) show profiles at constant CM time in the compression stage (for (a) 
the time steps are 0,20,40, 100, for (b) 0,200,400, 1000), (c,d) the subsequent expan- 
sion (for (c) the time steps are 120, 140, 240, for (d) 1200, 1400, 2400). For the sake of 
clarity, profiles are alternatingly shown as full and dotted lines. 

Fig. 10: The same as in Fig. 9 for the SHASTA (A = 0.4) with (a,c) Ax = 0.025 R 
and (b,d) Ax = 0.0025 R. 
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Fig. 11: The effect of reducing the antidiffusion fluxes in the SHASTA, (a-c) as in Figs. 5 
(d-f), (d-f) as in Figs. 6 (d-f). 

Fig. 12: The effect of using the physical velocity of sound in the signal velocity estimates in 
the relativistic HLLE (c,d) as compared to the correct solution (a,b) (Figs. (a,b) are identical 
to Figs. 5, 6 (b)). 
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